% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Function for loading bond data (input N in months)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = loadBondData(N)

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOAD GSW DATA
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NN = 360;
if (N>NN), error('Max maturity is 30 years'); end

tmpNom = readtable('../input/feds200628.csv', 'Range', 'A10');

idx = ~isnan(tmpNom.BETA0(:,1));
dates = datenum(tmpNom.Date);
dates = dates(idx);
Beta = [tmpNom.BETA0, tmpNom.BETA1, tmpNom.BETA2, tmpNom.BETA3];
Beta = Beta(idx,:);
Tau = [tmpNom.TAU1, tmpNom.TAU2];
Tau = Tau(idx,:);
mdy = tmpNom.Date;
mdy = mdy(idx);

n = (1:1:NN)/12;
novertau1 = (1./Tau(:,1))*n;
novertau2 = (1./Tau(:,2))*n;

yNom1 = Beta(:,1);
yNom2 = Beta(:,2).*(1-exp(-novertau1))./novertau1;
yNom3 = Beta(:,3).*((1-exp(-novertau1))./(novertau1)-exp(-novertau1));
yNom4 = Beta(:,4).*((1-exp(-novertau2))./(novertau2)-exp(-novertau2));
yNom = yNom1+yNom2+yNom3+yNom4;

pNom = -(ones(numel(mdy),1)*n*12).*yNom/1200;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Monthly Data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
breakM = find(abs(diff(month(mdy))))+1;
M = numel(breakM);
breakM = [1; breakM; numel(mdy)+1];
%
gsw.m.dates = NaN(M+1,1);
gsw.m.prc = NaN(M+1,N);
gsw.m.yld = NaN(M+1,N);
%
for i = 1:M+1
    tmpP = pNom(breakM(i):breakM(i+1)-1,:);
    tmpY = yNom(breakM(i):breakM(i+1)-1,:);
    gsw.m.dates(i) = 100*year(mdy(breakM(i)))+month(mdy(breakM(i)));    
    gsw.m.prc(i,:) = tmpP(end,1:N);
    gsw.m.yld(i,:) = tmpY(end,1:N);
end
%
gsw.m.period = 12;
gsw.m.mats = (1:size(gsw.m.prc,2))/gsw.m.period;
%    
gsw.m.fwd = gsw.m.yld;
gsw.m.fwd(:, 2:end) = 100 * gsw.m.period * (gsw.m.prc(:, 1:end-1) - gsw.m.prc(:, 2:end));
%    
gsw.m.ret                = NaN(size(gsw.m.yld));
gsw.m.xret               = NaN(size(gsw.m.yld));
gsw.m.dret               = NaN(size(gsw.m.yld));
%    
gsw.m.ret(2:end, 2:end)  = 100 * gsw.m.period * (gsw.m.prc(2:end, 1:end-1) - gsw.m.prc(1:end-1, 2:end));
gsw.m.ret(2:end,1)       = gsw.m.yld(1:end-1, 1);        
%
gsw.m.xret(2:end, 2:end) = gsw.m.ret(2:end, 2:end) - gsw.m.ret(2:end, 1) * ones(1, size(gsw.m.ret(:, 2:end), 2));
gsw.m.dret(:, 2:end)     = gsw.m.ret(:, 2:end) - gsw.m.ret(:, 1:end-1);


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QUARTERLY DATA
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
breakM = find(abs(diff(quarter(mdy))))+1;
M = numel(breakM);
breakM = [1; breakM; numel(mdy)+1];    
%
gsw.q.dates = NaN(M+1,1);
gsw.q.prc = NaN(M+1,N);
gsw.q.yld = NaN(M+1,N);
%
for i = 1:M+1
    tmpP = pNom(breakM(i):breakM(i+1)-1,:);
    tmpY = yNom(breakM(i):breakM(i+1)-1,:);
    gsw.q.dates(i) = 100*year(mdy(breakM(i)))+quarter(mdy(breakM(i)));    
    gsw.q.prc(i,:) = tmpP(end,3:3:3*N);
    gsw.q.yld(i,:) = tmpY(end,3:3:3*N);
end
%
gsw.q.period = 4;
gsw.q.mats = (1:size(gsw.q.prc, 2))/gsw.q.period;

gsw.q.yld = -100 * gsw.q.prc./(ones(size(gsw.q.prc,1),1) * gsw.q.mats);

gsw.q.fwd = gsw.q.yld;
gsw.q.fwd(:, 2:end) = 100 * gsw.q.period * (gsw.q.prc(:, 1:end-1) - gsw.q.prc(:, 2:end));

gsw.q.ret                = NaN(size(gsw.q.yld));
gsw.q.xret               = NaN(size(gsw.q.yld));
gsw.q.dret               = NaN(size(gsw.q.yld));

gsw.q.ret(2:end, 2:end)  = 100 * gsw.q.period * (gsw.q.prc(2:end, 1:end-1) - gsw.q.prc(1:end-1, 2:end));
gsw.q.ret(2:end,1)       = gsw.q.yld(1:end-1, 1);        

gsw.q.xret(2:end, 2:end) = gsw.q.ret(2:end, 2:end) - gsw.q.ret(2:end, 1) * ones(1, size(gsw.q.ret(:, 2:end), 2));
gsw.q.dret(:, 2:end)     = gsw.q.ret(:, 2:end) - gsw.q.ret(:, 1:end-1);

data.gsw = gsw;

end